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Numerical solutions of the Landau-Lifshitz-Gilbert micromagnetic model incorporating thermal 
04 ' fluctuations and dipole-dipole interactions (calculated by the Fast Multipole Method) are presented 

for systems composed of nanoscale iron pillars of dimension 9nm x 9nm x 150 nm. Hysteresis loops 
generated under sinusoidally varying fields are obtained, while the coercive field is estimated to be 
1979 ± 14 Oe using linear field sweeps at T — OK. Thermal effects are essential to the relaxation 
C*") , of magnetization trapped in a metastable orientation, such as happens after a rapid reversal of an 

external magnetic field less than the coercive value. The distribution of switching times is compared 
to a simple analytic theory that describes reversal with nucleation at the ends of the nanomagnets. 
Results are also presented for arrays of nanomagnets oriented perpendicular to a flat substrate. Even 
at a separation of 300 nm, where the field from neighboring pillars is only ~ 1 Oe, the interactions 
have a significant effect on the switching of the magnets. 



I. INTRODUCTION 



Several emerging technologies will incorporate fabricated magnets that are small enough to contain only a single 
magnetic domain. These nanoscale magnets will be essential for smaller components, lower power consumption, and 
completely new applications in fields such as information storage, integrated circuits, sensor technology, and micro- 
electromechanical systems. Successful implementation of these technologies requires a fundamental understanding of 
O , the dynamics of the internal magnetic structure of the nanoscale magnets on time scales ranging from the nanoseconds 
associated with gigahertz applications to the years over which magnetic information storage must be stable. 
| An essential factor in many of these applications, especially in information storage, is the free-energy barrier that 
. separates two antiparallel orientations of the magnetization. This free-energy barrier can be surmounted using thermal 
energy momentarily "borrowed" from the surroundings, and often device engineers strive for barriers of at least 40fceT 
to make this thermal bit-flipping a rare event. Here T is the absolute temperature and /cb is Boltzman's constant. 
The particle can be magnetizes in a specified direction by applying a field strong enough to remove the free-energy 
barrier. The smallest field sufficient to do this at zero temperature is called the coercive field H c . Fields smaller than 
the coercive field cannot change the magnetization orientation deterministicalhz. but they do cause lower free-energy 
barriers and make thermal crossing more probable. In fact, hybrid recording^ uses lower-than-coercive fields and 
. thermal barrier crossing to write data in high-coercivity magnetic media. 

Two examples of particular thermally activated barrier crossings are shown in Fig. [l| Red is associated with the 
magnetization orientation antiparallel to the applied field, which is metastable. Blue is associated with the parallel 
magnetization, which is thermodynamically stable. Green and yellow are intermediate between the two orientations. 
For times before those shown in the figure, the magnet stays in the metastable state with end caps of non-uniform 
magnetization induced by pole avoidance. Eventually thermal fluctuations in the end caps carry the system over the 
free-energy barrier. After the barrier crossing, the entire magnet changes quickly to the stable equilibrium orientation. 

Often, nanoscale magnets are assumed to be uniformly magnetized bodies. For instance, Garcia-Palacios and 
LazaroB considered the stochastic trajectories of isolated magnetic moments using the Landau-Lifshitz-Gilbert equa- 
tion subject only to applied fields and a uniaxial anisotropy, which supplied the free-energy barrier. They measured 
the susceptibility of the magnetization in small sinusoidal probe fields. Multiple crossing of the free-energy bar- 
rier, a frequently noted consequence of the gyromagnetic motion in single-spin models, was also observed for their 
uniform-magnetization model. 

When at least one dimension of the nanomagnet is greater than the exchange length of the material, nonuniform 
reversal modes become energetically possible. To consider such nonuniform magnetization reversal we use micromag- 
netic simulation with a large number of points inside the magnet. Specifically, nanomagnets with an aspect ratio 
of approximately 17 are considered throughout this work, and arc referred to as pillars. The numerical approach 
is discussed in Sec. |fij with further details found in the appendices. Results for the estimation of the coercive field 



and the statistics of magnetization switching are discussed in Sec. III. A simple model for nonuniform switching is 
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discussed in Sec. |y|. Section |y| considers the effect of interactions between pillars in a one-dimensional array. A brief 
summary of the results is given in Sec. VI. Some preliminary results of this study were presented in Refs. |3|-|5j. 

II. NUMERICAL DETAILS 

The basic approach of micromagnetic modeling is to consider a system of coarse-grained magnetization vectors 
M(?i), with ?i indicating the location in space of the z-th spin. The vectors are assumed to have a fixed magnitude 
M s , corresponding to the bulk saturation magnetization density of the material. This is a valid approximation for 
temperatures well below the Curie temperature.Q The evolution of each magnetization vector is governed by the 
damped precessional motion given by the Landau-Lifshitz-Gilbert (LLG) equatioroJ 

dM(Fj) _ 70 a - 



di 



r ^M(r l ) x (fift) - ]j£M(?0 x H(?o) (1) 



where the electron gyromagnetic ratio is 70 = 1.76 x 10 7 Hz/Oe,u and a is a phegomenological damping parameter. 
For the sign of the undamped-precession term we follow the convention of Brown.0 Here a tilde is used to distinguish 
dimensional quantities from their dimensionless counterparts, which will be introduced later. 

The dynamics are controlled by H(?j), the local field at the i-th position, which, in general, is different at each 
lattice site. The local field mediates all of the interactions in the system with the contributions combined via linear 
superposition, 

H(?j) = H z (fO + H c (r,) + H d (f,) + H a (?i) + H n (?i) , (2) 

where H z (?i) represents the externally applied field (Zeeman term), H e (Fj) is due to exchange effects, Hd(?j)is the 
dipole field, H a (?i) is due to crystalline anisotropy, and H n (?i) is the random field induced by thermal noises The 
estimates used for each of these fields are given below, except for H a (fj), which we have taken to be zero. In this 
article, a ~ emphasizes that the symbol represents a quantity with units, the corresponding dimensionless quantity 
is written without the tilde. Material parameters, such as M s , always have units when appropriate. 

We have chosen material parameters to match those of bulk iron; a summary of all parameters used here is given 
in Table 1. The saturation magnetization density isJkL = 1700 emu/cm 3 ,E_2l while the exchange length, the length 
over which M can change appreciably, is l c — 3.6 nm£lo The damping parameter Qf-Jaas proven difficult to measure 
or estimate from ab initio considerations, and may even depend on numeric details.113 Here we have chosen a = 0.1 
to represent the under-damped behavior usually assumed to exist in nanoscale magnets. In the other extreme, the 
over-damped limit, the gyromagnetic motion intrinsic, to the LLG equation is suppressed, and the system can be more 
efficiently simulated using Monte Carlo methods.oll3 

The magnetization was discretized on a cubic lattice with discretization length Af = 1.5 nm chosen to give sufficient 
resolution across the cross section of the pillar. The time discretization of At — 5 x 10~ 14 s was chosen to avoid 
numerical instabilities in the integration. The details of the integration scheme are given below. The equation of 
motion was cast into dimensionless quantities by considering the rescaled length r = r// 0f _±ime t — 7oAf s i, field 
H = H/M s , and magnetization density M = M/M s . This rescaling is facilitated in cgs units, EJ where magnetic fields 
can be normalized by the magnetization density without the inclusion of any constants. 

In a continuum model, exchange represents local differences in the alignment of-jthe magnetization. The contribution 
to the local field from exchange interactions can be approximated by / 2 V 2 M(f )El which has been implemented in the 
simulations by 

H e (r i )=^V l-6M(r i )+ £ M(r< + d) J , (3) 

where the summation is over the six nearest neighbors of Fj. Here the exchange length is defined in terms of the 
exchange energyliS E c = —(1^/2) j dr M • V 2 M. The dimensionless form of the exchange contribution to the local field 
is 

He ^)= (at) 2 f- 6M K>+ E M(r, + d)l . (4) 



|d|=Ar 
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The computationally most intensive part of determining the local field is finding the contribution from dipole-dipole 
interactions, Hd(r). Efficient calculation, which uses the Fast Multipole Method (FMM),E3 is quite involved and is 
discussed in Appendix 1. 

Thermal fluctuations are included in the LLG equation of motion by inclusion in the superposition of a random 
field H n (r) with Gaussian fluctuations whose first moments are zero and whose second moments obey the fluctuation- 
dissipation relatione! 

(H^ihMn,' = (t-t>) 8^8 i>v , (5) 

where H nfi indicates one of the Cartesian components of H n . Here V — (Af) 3 is the discretization volume of the 
numerical implementation and 8^^ is the Kxonecker delta representing the orthogonality of Cartesian coordinates- 
This result was derived for isolated particlesfl, and interactions between the discretization volumes can be important.tS 
To make the simulations tractable, the effect of interactions on the thermal noise has been neglected. Here 8(t — t') is 
the Dirac delta function, and its dimensions [8(t — t')] = 1/s are important. After the transformation to dimensionless 
quantities, the fluctuation-dissipation result is 

(H^Ti^H^ir'^t')) = e8{t-t')8^,8i, v , (6) 

with the dimensionless strength of the stochastic field given by 

2ak B T 



(7) 



where the untransformed discretization volume and saturation magnetization are both used in cgs units to yield the 
dimensionless result. The dimensionless form of the LLG equation is 

= Y^mr*) x [H(r,) - oM( ri ) X H(r,)] . (8) 

This dimensionless stochastic differential equation is used for the numerical integration. 

Since the LLG equation conserves the magnitude of the magnetization density, each integration step amounts to a 
rotation. The size of the tangential displacement at each integration step is given by 

R(r<, t) = ——M(Ti,t) x [(I(r 4 , t) - aM(r„t) x I(r 2 , *)] , (9) 
1 + or 

where I is the "impulse" over the integration step defined below. The magnetization after the integration step is given 
by 

M(ri,t) +R(r 4 ,i) , . 

M(f "' + A "= Vl + R^U • (10) 

which ensures the conservation of unit magnitude for M. If only the deterministic dynamics were important a high- 
order integration method could be used, but to correctly include the thermal noise a lower-order method must be 
used. Using first-order Euler integration, the impulse including the fluctuating field is 



I(r,,t) = [Hafa.tJ+Hdfa.tJ+Hefa.t)] At +V7Atg(r,,t) , (11) 

where g(r, t) is a random vector with each component chosen independently from a Gaussian distribution of zero mean 
and unit variance. This result, including the %/ At contribution characteristic of integration of stochastic processes, is 
explained in Appendix 2. 



III. ISOLATED NANOMAGNETS 



The numerical micromagnetic simulation methods described in the previous section have been applied to model 
nanoscale magnets inspired by those recently fabricated using a combination of chemical vapor deposition and scanning 
tunneling microscopy techniques by Wirth, et aloa The technique has been used to produce arrays of nanoscale 
magnetic particles with diameters down to about 10 nm and lengths from 50 to 250 nm. The reversal of pillars with 
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diameters more than twice the exchange length, such as these, have been_£pund numerically to proceed by a mode 
where the magnetization is not constant across the diameter of the pillar .EETlj The model particles considered here 
are rectangular prisms 9 nm x 9 nm x 150 nm, which require N = 4949 sites on the computational lattice. The long 
axis is taken as the z-axis, and the cross-sectional area is chosen comparable to that for a 10 nm diameter pillar. The 
external field is always applied along the z-axis - the easy axis induced by shape anisotropy. As mentioned in Sec. [fi|, 
the material parameters were chosen to ma tch .those of bulk iron. These are consistent with measurements on the 
experimentally produced iron nanoparticles .113113 



A. Coercive Field 

Images of the z-component of the magnetization, M z , for such a particle as described above are shown in Fig. ||(a) 
for several fields during a hysteresis- loop simulation (in a one-quarter cut-away view). The color is a linear scale of 
M z , shown in the legend of Fig. [j] with the most negative values at the bottom and the most positive values at the top. 
For hysteresis-loop simulations, the field is always oriented along the z axis. Starting from its most positive extreme, 
the value is varied sinusoidally. For the simulation shown in Fig. || a period of one nanosecond, a field amplitude of 
4000 Oe, and a temperature of K were selected. The response of the initially upward-directed magnetization after 
the applied field is oriented downward can be clearly seen. First, end caps associated with pole avoidance form at 
both ends. These end caps are regions of high curl in the magnetization density, C(r) = V x M(r). The end caps are 
obvious in Fig. 0(b), where sign (C z (r, t)) |C (r, t) | is shown on a linear color scale, so that right (left)-handed end 
caps are colored red (blue). At zero temperature, the two end caps then grow symmetrically until they meet near 
the midpoint of the particle. While they grow, the area of large curl is concentrated at the interface between the 
volume where the magnetization has already aligned parallel with the applied field and the volume where it remains 
antiparallel to the field. Each interface is composed of two regions with opposite z-component of the curl, with a 
curling type of singularity in the center of the pillar. These two regions are separated by a layer where the z-component 
of the curl is close to zero, and the magnetization vectors in the xy-plane have a negative divergence located at the 
center of the pillar. Because the end caps have opposite curls, in the region where the end caps come into contact the 
curl changes abruptly from large positive to large negative values. Some time is required for this defect to disappear. 
This reversal of pillars by nucleation of reversed volumes at the ends followed, by growth of the reversed regions is 
consistent with minimization of the micromagnetic energy at zero temperatureBEj and Monte Carlo simulations at 
finite temper ature,tfl both under quasistatic field-sweep conditions. 

The hysteresis loop associated with this simulation is shown in Fig. ||, with the field values shown in Fig. || indicated 
by large tick marks. Two periods are shown for the Ins loop. Aside from differences due to the initial alignment 
of the magnetization and the presence of the defect after a complete loop, the reproducibility at zero temperature 
is excellent. Hysteresis loops under the same conditions, but for periods of 2 and 4 ns, are also shown. These 
simulations show that the hysteresis loop becomes more square for longer periods, indicating that the magnetization 
is not following the applied field in a quasi-static fashion. This is quite reasonable given that the simulated loops 
correspond to microwave frequencies, but it makes the infinite-period coercive field difficult to measure. 

The average energy density, along with its contributions, as functions of time appear in Fig. || for the 4 ns hysteresis 
loop. The energy density for zero crystalline anisotropy is calculated as 

1^/1 \ 

E = -jj E M ( r <) • o [Ho(r ° + Hd(ri)] + Hz(r ° ' (12) 

»=i ^ ' 

where N is the number of lattice sites. As the applied field begins to decrease and then become negative, it is the 
Zeeman energy that changes the most while the exchange and dipole-dipole energies remain nearly constant. Near 
the coercive field, the exchange energy rapidly increases as regions of reversed magnetization, and the interfaces 
associated with them, form at the ends of the pillar. Once these reversed regions become large enough, they grow 
spontaneously and the Zeeman energy rapidly decreases, while the exchange energy from the interface between reversed 
and unreversed regions remains roughly constant. The growth of the reversed regions does not take long, and the 
exchange energy quickly dissipates when the regions from the two ends merge at the middle and the reversal is 
complete. 

The information in Fig. |] can be used for extracting the coercive field of the nanomagnet. For infinitely slow 
variations in the applied field, dE/dt is proportional to dE / dH a using the chain rule. Thus portions of the figure with 
positive slope of the total energy correspond to fields for which the system is trapped behind an energy barrier, while 
portions with negative slope correspond to fields for which the system responds deterministically by finding a new 
energy-minimizing magnetization configuration. Fields for which the energy has a maximum should then correspond 
to the coercive field. Dynamic effects cause the field at which this maximum occurs to depend on the frequency. 
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To estimate the static coercive field we find H C (R) for fields that vary linearly with time, H(t) = —Rt for t > 0, 
and then extrapolate to the H C (R = 0) value. The energy density is presented as a function of the absolute value of 
the applied field in Fig. [5] for rates ranging from R = 250Oe/ns to R — 10000 Oe/ns. The coercive fields estimated 
from the energy maximum are shown in the inset where the error bars are estimated from the second derivative near 
the energy maximum. The coercive field clearly decreases with R. Extrapolation of a weighted least-squares fit yields 
a coercive field of H c = 1979 ± 14 Oe at R = 0. This is a more accurate estimate of the static coercive field than those 
taken directly from hysteresis-loop simulations. 

B. Switching in Constant Field 

At finite temperature the magnetization can reverse, even when the applied field is weaker than the coercive field. 
In this case, thermal fluctuations carry the magnetization past the free-energy barrier and on to the new equilibrium 
configuration. Two examples of such thermally induced switching are shown in Fig. [l] for the conditions Hq = 1800 Oe 
and T — 20 K. For small positive times (not shown) the volumes associated with the end caps fluctuate due to the 
thermal noise. The reversal appears to proceed by a nucleation process, with the end caps serving as the seeds for 
heterogeneous nucleation. The initiation of the switching from the end caps is similar to results for reversal induced 
by applied fields greater than H c in experiments,E3 and simulations with T=0 K where reversal begins at the end 
cap soo or corners E3 

In classical droplet theory, in which boundaries such as the sides of the magnet are ignored, nucleation of the 
equilibrium magnetization is governed by the competition between the favorable Zeeman energy due to alignment 
with the applied field and the unfavorable exchange energy due to the interface with the majority volume of the 
misaligned orientation. The exchange energy dominates for small droplets, which tend to shrink. Large droplets tend 
to grow because the Zeeman energy dominates. These two regimes are separated by a critical droplet size, the saddle 
point associated with the free-energy barrier, where the tendencies toward growth and shrinkage are balanced. In 
nanoscale magnets the situation is more complicated because the free energy cannot be easily separated into surface 
and volume contributions. Nevertheless, the important aspect remains that a free-energy barrier exists that must 
be crossed for reversal to begin. The end caps fluctuate until a succession of highly unlikely fluctuations carries the 
end-cap configuration past the free-energy barrier. After that, growth of the reversed region is energetically favorable 
and occurs rapidly, on the order of 0.2 to 0.3 ns for the pillars. 

For long pillars, nucleation events at the two ends occur approximately independently of each other. In addition, 
the growth of a supercritical region takes a significant amount of time: enough time for the other end cap to have 
a reasonable probability to nucleate and begin to grow. One example of both ends nucleating at different times is 
shown in Fig. [l](a), while Fig. |l|(b) shows nucleation of one end that grows to switch the magnetization of the entire 
pillar. For the present applied field the latter situation is somewhat rare, from the theory presented below it can be 
inferred that it occurs in about 10% of simulated switches. From Fig. |l|(b) it can also be seen that the nucleation of 
the two end caps is not completely independent. As one reversed region grows along the pillar, the free-energy barrier 
at the other end cap decreases. Eventually the free-energy barrier for the second end becomes zero and it will also 
begin to grow. If a less negative field is applied, the free-energy wells of the metastable states will be deeper and the 
first reversed region will have to grow further along the pillar before this occurs. 

The mean switching time, F sw , defined here as the first time when M z = 0, is shown in Fig. ^(a) as a function of 
applied field for T = 20 K. The difference in the error bars reflects the amount of statistical sampling at the different 
conditions. The dependence of t sw on temperature is shown in Fig. ^(b) for applied fields of Ho — 1850 Oe; fields 
are always applied along the pillar. While there is a clear dependence on field and temperature, we find no clear 
functional form for the dependence. Notably, the exponential dependence of t sw on inverse temperature, expected 
when the barrier is high, is not seen. This implies that the free-energy barrier is low. In fact, this picture is consistent 
with the switching statistics described next, where the dynamics are consistent with a biased walk in a shallow well. 

The statistics of magnetization switching are measured by P no t(t), the probability of not switching before time t. 
The simulation results for P no t{t) for 200 switches at Ho = 1800 Oe and T = 20 K are presented in Fig. ^(a). The 
result is clearly not exponential, the expected functional form when the free-energy barrier is large. Experimental 
results for single-domain magnets, in which P no t(t) is not exponential, have been reported recentlyES 

A simple theory that assumes completely independent nucleation at the two ends of the pillar produces a reasonable 
description of the simulation result. The theory for P no t(t) requires two parameters. The first is the constant nucleation 
rate, X, for the formation of a supercritical reversed region at one isolated end of a pillar. The second is the rate, v, 
at which a single growing region changes the normalized magnetization of the pillar. Labeling the two independent 
nucleation times at the top and bottom i t and tb, respectively, the positive values of the magnetization of the pillar 
can be described by 
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1 t<tl 

M z (t) = { l-v{t-h) h<t<t 2 , (13) 

1 - v(t - 1±) - v(t - t 2 ) t 2 <t 

where t\ = min(t t ,tb) and t 2 — max(i t , t^). We define switching to occur at M z (t svr ) = 0, when roughly one-half the 
volume of the nanomagnet is oriented in the equilibrium direction. The switching can occur either through one or 
both ends nucleating, and the actual switching time for a particular (tt)ib) is given by 

W = minf2t +*i J *o + ^Y^J , (14) 

where t = 1/(2?;) is the earliest time at which the pillar can switch. For times t < to the probability of switching is 
zero, even if both ends nucleate immediately after the reversal. The situation in the t^-tt plane is shown in Fig. ||, 
where the dashed curves correspond to two different conditions of constant t sw , t sw /t =3/2 and 3. The solid lines divide 
the plane into regions with different switching histories. In the triangle near the origin (I), switching occurs for times 
to < t sw < 2io, which can happen only if both ends nucleate. Between the parallel solid lines (II) switching occurs 
by double nucleation, while outside them (III & IV) a single nucleation causes the switching. Since each nucleation 
process has a constant rate, the corresponding probability density function for the nucleation time is exponential, 
Xexp {—It). Integrating along curves of constant t sw , the probability of not switching before time t is found to beB 

( 1 t<t 
Pnot (t) = { e- 21 ^-'") [1 + 21(t - 1 )} t <t<2t . (15) 

[ e -2X(t-*o) [j + 2Jt ] 2t < t 

We will refer to this as the "two-exponential" decay model. In the limit of infinitely fast growth, to — 0, only regions 
III and IV remain, and this becomes a simple exponential decay associated with thermally activated crossing of a 
free-energy barrier. For to > 0, P no t(t) is quadratic in t — to for t — > ij, has a discontinuity in the slope at t = 2to, 
and is exponential for t > 2to- 

The fit of the two-exponential decay model of Eq. (|l5| ) to the simulation results appears in Fig. 0(a) as the dashed 
curve. The dotted curve is a fit to the error function 

PcAt) = \ - ±Erf (l(t - t )) , (16) 

which was chosen as an empirical form because it also has two adjustable parameters, chosen here to have a similar 
interpretation as the theory given above. The parameters for both theoretical forms were determined by matching 
the first and second moments of the theoretical form to the moments of the simulation data, and the parameters are 
2 ps 4.049, t » 0.527 and 2 ps 4.176, t ps 0.790. For these conditions, both forms, Eqs. © and ©, seem to 
fit the simulation data about equally well. The disagreement with Eq. ( |l5| ) probably occurs because the condition 
that the free-energy barrier is much greater than k^T is not met. When the free-energy barrier is much less than 
k^T, the magnetization should be more like a biased random walk. In that case the most noticeable effect of thermal 
fluctuations would be a symmetric smearing of the switching times around t sw . 

The question of statistical sampling is also important. Specifically, one expects to have undersampling of the 
population because the finite-sampling P n ot(0 goes to zero faster than the infinite-sampling P n ot(i) with increasing 
t. In fact, for times less than the maximum observed switching time in the long-time tail of the distribution, the 
error function underestimates the population, while the two-exponential theory overestimates it. To this extent, it 
is possible that data with better sampling may shift towards better agreement with Eq. (^). Indeed, experimental 
switching probabilitiesEa in single-domain magnets have been observed to have exponential tails, which is inconsistent 
with Eq. ©. 

Better statistics and a higher free-energy barrier are needed to validate the theoretical P no t(t) given in Eq. (jjl]). 
For the present model, better statistics would be prohibitive; the 200 switches considered here took approximately 20 
weeks to generate using 4 processors on an Origin 2000. A less realistic but numerically faster simulation, in which 
the pillars are modeled as one-dimensional stacks of cubes, is presented in the next section. The advantage of this 
simplified model is that it allows many more switching events to be observed in a reasonable amount of computer 
time, even with much higher free-energy barriers. This improves the statistics for the very rare events with i sw 3> 2to 
and t sw — > t , which are important to distinguish between alternative theoretical forms for P n ot(*)- 
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IV. A SIMPLE MODEL 



In order to generate a sufficiently large sample of switching events in nanoscale magnets, we consider a simpler 
model in which, the pillar consists of a one-dimensional stack of cubic cells, which is the model introduced by Boerner 
and Bertram.Ej The cubes have a linear size Af = 2£ c , here 5.2 nm, since the material parameters for iron are used 
(Rcf. used parameters for nickel). To keep the aspect ratio of the nanopillars as close as possible to those considered 
above, a pillac, composed of 17 cubes is used. For this simple model, the local field due to the dipole interactions is 



calculated as 



H' d (f;) = (Ar) 3 £ —A_ _p! — , (17) 

where Fy is the displacement vector from the center of cube i to the center of cube j, and is the corresponding 
unit vector. The factor of volume (Af) 3 results from integrating over the constant magnetization density in each cell. 
This factor can be combined with rfj, so that the denominator depends only on Xy = ry/(Af), which is a vector 
of integers denoting the difference in cell indices. It can readily be seen that for a uniformly magnetized pillar the 
summation over index integers leads to different values of H[j(r) for different choices of discretization length along 
the long axis. This happens because self-contributions have been ignored and the far-field result has been used for 
neighboring cells. However, since this model is being used to investigate only the statistics of switching, and not for 
estimating physical values, the results should be qualitatively correct. Using the scaling to dimensionless units, the 
approximate dipole field is 

H^)^E 3f " (? " ,M( 3 3)) " M(rj) - (18) 

We have verified that our implementation of the model agrees with that of Boerner and Bertram by reproducing the 
observed coercive field for the nickel pillars discussed in Ref. Because of the different approximations, the coercive 
field for the present iron pillars is about 1500 Oe. 

The probability of not switching for 2000 switches for this model of iron pillars is shown in Fig. Wb), along with the 
two theoretical forms, Eqs. ([l5]) and (|l6|), which have been fit using the same procedure as in Sec.jnJ. The kink in the 
two-exponential theoretical form at 2to is quite noticeable, but cannot be seen in the simulation data. The kink has its 
origin in the complete suppression of nucleation for negatives times, and the absence of the kink in the simulation data 
may stem from the finite chance of nucleation at negative times due to the way the field is reoriented, as described in 



Sec. [II A. Another possibility is that the interactions between the ends smooth out the difference between the one- 
and two-nucleation decay modes. Despite the kink, the theoretical form Eq. ([l5]) does a good job of describing both 
the exponential tail and the rounding at early times of P n ot{t)- It is clearly superior to the error- function form, as 
well as to a displaced, single exponential (not shown). 



V. ARRAYS OF NANOMAGNETS 



The long-ranged dipole-dipole interactions that contribute to the shape anisotropy of nanoscale magnets also cause 
interactions between nanomagnets. This is especially true in most potential applications, where miniaturization will 
drive devices to high densities. In addition, the interactions between nanomagnets in arrays could be the basis of 
device applications, prototypes of which have already been investigated. EJEiJ Regular arrajrs .pf. nanomagnets have 
already been used experimentally to provide magnetic signals strong enough to be measuredtllEj in the experiments 
our nanomagnets are modeled after. Our simulations show that even for very wide spacings the magnetic interactions 
between nanomagnets have significant effects on the switching properties. 

To allow for analytic treatment, we consider only the two leading contributions to the local dipole field Hd(r) 
from the dipole and quadrupolc moments of the source magnet. The former contributes uniformly throughout the 
observation volume, while the latter changes linearly in each Cartesian direction. Specifying the distance between 
pillars as d and considering the dipolar and quadrupolar moments of the source magnet, and M®, respectively, 
the leading contributions to the observed demagnetizing field are H^(r) = -e z M°/d 3 and H^ 2) (r) = -9(2ze z - 
xe x — ye y )M2 /(4d 5 ), respectively. Simple expressions for the multipole moments can be calculated in the following 
way. Assume a pillar with regions of uniform magnetization oriented in the -l-e direction in the middle and in the 

e direction at the ends, as shown in Fig. |9|. With the total pillar length L, the top-region length £ t , and the 
bottom-region length £^ all measured along the long axis of the pillar, the dipole moment for a reversing pillar is 
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and the quadrupole moment is 



= M S A (L - 2£ t - 24) (19) 



M§=2M s A(%-%-L(i t -e h )) , (20) 



where A is the cross-sectional area of the pillar. The dimensionless form is achieved by setting Ms to unity and using 
dimensionless lengths. 

Using the simplified expression for the dipole moment in Eq. ([l9|), the contribution from one upward-magnetized 
pillar at the nearest-neighbor position in the array is • e z w — 1 Oe. The quadrupole contribution from Eq. ( p0| ) at 

(2) 

the end of one nanomagnct that is the nearest neighbor to a switching pillar with £ t = L/2 is • e z ss 0.4 Oe. This 
is a rough estimate of the maximum interaction through the quadrupole moment, and the sign indicates a tendency 
for neighboring pillars in the array to switch at opposite ends. 

As a simple initial investigation, we consider a linear array of four of the rectangular nanomagnets described in 



Sec. [II, with the nanomagnets oriented perpendicular to the substrate. Since dipole-dipole interactions within each 
nanomagnet are calculated using the Fast Multipole Method described in Appendix 1, the multipole moments for 
each nanomagnet are readily available. These moments can be used to quickly calculate the interactions between 
nanomagnets in the array, under the constraint that the far-field description is appropriate. To ensure this we have 
considered only a spacing between pillars of two pillar lengths, 300 nm. Note that this situation is quite different from 
Fast Fourier Transform approaches, in which the calculation must be carried out on a lattice that also fills the space 
between the nanomagnets£3 To do that practically, the space between the magnets must be kept small. 

To study arrays of nanomagnets, systems consisting of four 9 nm x 9 nm x 150 nm parallel pillars arranged in 
a line perpendicular to their long axes and spaced 300 nm apart were simulated using the Fast Multipole Method 
truncated at p = 3. Hysteresis loops with periods on the order of a few nanoseconds for the individual pillars in the 
array look similar to those for isolated pillars in Fig. [| with no observable difference between pillars on the outside 
and those on the inside of the array. The symmetry equivalence for the two pillars on the inside, and for the two 
pillars on the outside, will be used throughout to double the statistical sampling. 

The probability of not switching for H = 1800 Oe and T = 20 K is shown for 40 array switches in Fig. |l0|(a) for 
pillars on the inside and outside of the array, as well as isolated pillars. No significant difference can be seen between 
the three curves. However, the coupling between the pillars can be seen in the difference between the P no t(t) for inside 
pillars with one or both nearest neighbor pillars switched, shown in Fig. |l^(b). Here, t is the time difference between 
t sv , and the last time a neighboring pillar switched. From this data it can be seen that of pillars with two neighbors, 
those with only one neighbor switched tend to switch earlier than these with both pillars switched. The effect is even 
more pronounced in simulations of the simple model at H^lOOOOeB 



VI. SUMMARY 



Numerical simulation of the Landau-Lifshitz-Gilbert micromagnetic model has been used to investigate spatially 
non-uniform magnetization switching in nanoscale magnets at the nanosecond time scale. We have focused on.ir.on. 
pillars 9 nm x 9 nm x 150 nm because such pillars and arrays have been constructed and measured experimentally£j : E£l 
The zero-temperature static coercive field has been estimated numerically to be H c — 1979 ± 14 Oe by finding the field 
of maximum energy for fields swept at constant rate, and then extrapolating to find the zero-rate estimate. Simulations 
of thermally activated magnetization switching are possible at fields well below the coercive value. For the pillars 
studied here, reversal occurs through nucleation at the ends of the pillars. The probability of not switching, P no t(t), is 
not well described by a delayed exponential, and a theoretical form, Eq. ( |l5|) , based on independent nucleation at the 
ends of the magnet is developed here. The agreement with results for intensive, fully three-dimensional simulations 
with an applied field near the coercive value are reasonable, but an ad hoc error function gives similar agreement. The 
agreement with Eq. Jl5| ) is much better when the field is well below the coercive value and the statistics are better, 
which currently we have only studied with less-intensive simulations that discretize the nanomagnet only along its 
long axis. 

The fully three-dimensional micromagnetics program has been developed for massively parallel computers and 
implemented using the Fast Multipole Method. These two features make it feasible to simulate nanomagnets in 
widely spaced arrays. Even though the interactions between the magnets are quite weak for the specific linear array 
considered here, there are significant effects on the statistics of the magnetization switching. Specifically, there is a 
dependence of P no t (t) for a given pillar on the orientation of the magnetization of its nearest neighbors in the array. 
The nature of the cooperative reversal mode observed here is a topic for future research. 
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APPENDIX 1 




Calculating the dipole-dipole interactions is the most intensive part of the numerical calculation. The magnetic 
potential approach used here involves defining a magnetic charge density pm(y) = — V • M(r). (We present this only 
in terms of the dimensionless quantities.) This charge was evaluated on a cubic lattice dual to that of M(r,), see 
Fig. [ll], using equally weighted two-point differences, specifically 

/<mUM> = - > \ | — signCh-tyMfo + h)-^ ) . (2:1) 

=\/3Ar/2 

where is a position on the dual lattice, is the sum over the corresponding corners of the cube from the direct 
lattice, and are the Cartesian unit vectors. If we were to apply Eq. ( |2~i| ) with just outside the magnetic material, 
it would give a nonzero charge there. To avoid this unphysical result, we have moved this charge to the surface by 
defining a surface magnetic charge density uyiiy) = s(r) • M(r), where s is the unit vector directed out of the surface. 
We have considered surface charges only on the surface of the model magnet, and we evaluate them at the centers of 
the squares defined by adjacent points on the surface of the direct lattice. The four corners are equally weighted so 
that 

a M (r d ') = j §'M(r d , +h), (22) 

|h|=V2Ar/2 

where runs over the four corners. The numerical approach of Eqs. ( pi"| ) and ( p2[ ) ensures that there is no net 
magnetic charge on the system as a whole. I— ■ 
The magnetic potential 0m (r) is found by integrating over both the volume and surface charges,E3 

Numerically, such an operation can be quite expensive since unsophisticated algorithms will require 0{N 2 ) operations, 
where Aj-is the number of lattice sites. An algorithm that remains reasonable for large systems is the Fast Multipole 
MethodO 

We have chosen to calculate the magnetic potential using the Fast Multipole Method (FMM) because it has several 
advantages over the more traditional Fast Fourier Transform (FFT) approach. The biggest difference between the two 
methods is that the FMM makes no assumptions about the underlying lattice, while the FFT method assumes a cubic 
lattice with periodic boundary conditions. One consequence of this assumption is that numerical models of systems 
without periodic boundary conditions require empty space around the magnet so that the boundary conditions do 
not affect the calculation. The FFT also requires the lattice to continue into regions of empty space that lie between 
elements of an array of magnets. By contrast, a FMM implementation only needs to consider volumes occupied by 
magnetic material. It does not need any padding. In addition to these advantages, the FMM is more efficient for 
large numbers of lattice sites. 

The popularity of the FFT approach stems from the fact that it takes O(NhiN) calculations to evaluate Eq. fl25|). 
The FMM has a larger overhead, but requires only O(N) operations to calculate the same potential.EJ This means 
an FFT approach to Eq. ( p3[ ) makes sense for small, cubic lattices, but that the FFM approach will be more efficient 
for large, irregular, or incomplete lattices. 

The FMM algorithm exploits the fact that 4>m at each lattice point can be expanded in terms of spherical harmonics, 

"u('-".~) = >: V ( I^r' + ^U , (24) 
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where the L\ terms can be used to represent the potential close to the lattice point, and the M? terms can be used 
to represent the potential far from it, but not both simultaneously. In this context, near means being closer to 
than any other lattice point, and far means distances more than twice the largest distance from the center of cell j to 
any of its boundary points. Following GreengardJiH we define the spherical harmonics 



Y*(6 t <p) = ^±^M P ^{ C0S 6)e 1 ^ , (25) 

with P? (x) the associated Legendre polynomial, and I = y— 1. Actually implementing this approach requires a 
truncation of the expansion in i at order p. We have found that the demagnetizing field for p — 3 is within 1% of the 
exact value for our simulations. 

Our implementation of the FMM algorithm starts by partitioning the model space into a system of cubes, an 
example is shown in Fig. [l2| The length of the side of the smallest cubes, which are separated by dotted lines, is 
the same as the discretization length. For our lattice system we have found that the most efficient choice is to use 
cubes centered on the direct lattice. For each cube, the multipole expansion coefficients of the far-field potential are 
calculated for the specific configuration of magnetic charges pm and cm 

Ml = [ dr [ PM + <tm] Q% j (0, if) , (26) 
Jv 

with the coordinates centered on the lattice site. Note here that g % is the distance from the center of the cell raised 
to the i-th power. For our cubic lattice, each quadrant of the cube is contributed by a different region of constant pm 
from the dual lattice (similarly for cm). With the geometry of the lattice fixed, the multipole expansion coefficients 
are easily calculated and summed to yield the total expansion coefficients for the lattice site. 

Each level of the hierarchy involves grouping cells into successively larger cubes that completely contain cubes at 
the lower level. The obvious hierarchy with each larger cube containing eight of the smaller cubes apparently works 
best. In Fig. [j^, cubes of the second level are separated by dashed lines and those of the third level by solid lines. 
Since the number of nodes in any direction along the simulation lattice is not restricted to a power of two, cubes do 
not always contain eight smaller cubes. The M\ for the larger-cube can be rapidly evaluated from those of the smaller 
cubes using the rule for translation of a multipole expansiorEH 

M i =EE or ^~'y : W ' , (27) 

k=0l=-k i 

where 0\ are the expansion coefficients for the smaller cube and the spherical coordinates (g,9,tp) here, are for the 
vector from the origin of the large cube to that of the small cube. The relations for the new factors aretll 

f-lV 

A\ = - [ ' (28) 

and 

jj = f ify<0 (29) 

1 \ 1 otherwise. 

The construction of the hierarchy terminates when the entire system is enclosed by a single cube. The partitioning 
can be generalized to noncubic rectangular prisms, but the restriction that the multipole expansion is only valid at 
distances larger than twice the diagonal of the rectangle will complicate the algorithm. 

A downward pass through the hierarchy involving three types of operations is required to construct 4>m as represented 
by the local expansion coefficients L\ for each cube. The first operation is the translation of the local expansion for 
the encompassing bigger cube (if one exists) that includes all the contributions ixpm its far field, i.e., those cubes on 
its level that are not its neighbors. This translation is accomplished by the ruleEZI 

ihEE O k / k ~^- A e "-%z{{0M , (30) 



k=i l=—k 



4. 



where 0\ are the local expansion coefficients of the larger cube, and the spherical coordinates (g, 9, if) here are for 
the vector from the origin of the smaller cube to that of the larger cube. A{ is defined in Eq. (ps|), whilcO 
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{-iy(-iy ifjfc<o 

Jij = { i£jk > Oand|fe| < |j| (31) 

(— l) 2 otherwise. 

The next operation incorporates the contributions from areas in the near region of the larger cube, but in the far 
region of the smaller cube. This is accomplished by transforming the multipole expansion of the source into a local 
expansion of the smaller cube usingO 

fc=0 l=-k A i+k U 

where the spherical coordinates are for the vector from the origin of the smaller cube to that of the larger cube andEzI 

jj,k = f if jk > o 

1 y (~1) 4 otherwise. 

The third type of operation is used for termination of the algorithm at the lowest level, where the near-region 
contributions for the smallest cubes must be evaluated exactly. Since a regular lattice is used, the point where (j>M 
is being calculated always lies at the corner of the neighboring cube of constant pm , and the contribution is simply 

(A02 W^VJWM, (34) 



V2J 2 

because only one quadrant of the region of constant /?m needs to be considered in this exact manner. The contribution 
for squares of surface charge that touch is 

Arsinh -1 (l)oM(r«i) ■ (35) 

Similar expressions can be calculated for general rectangular prisms. 

The FMM is an efficient way of determining the magnetic potential </>m(i"i) associated with a particular configuration 
of a dipole field. The local observed magnetic field due to the other dipoles is obtained using 

Hd(rO = -V<Mr*) (36) 

from potential theory. Numerically, this gradient was estimated with a centered difference of the nearest-neighbor 
sites, except on the boundaries, where forward or backward differencing was used. Note that the discretizations of all 
operators have been chosen for consistency between surface and volume charges such that the model magnet can be 
padded by volumes of lattice sites with M = without changing the results. 

The FMM was implemented using C++ and named Hierarchy. h. Fundamental to the implementation is the 
class sh_expansion whose instances are a (p + l) 2 array of complex numbers representing either local or multipole 
expansion coefficients. The class includes methods for indexing expansion coefficients within an expansion, evaluating 



the expansion at a specified point, and transforming an expansion using Eqs. (|27|), (30), and (p2[). In addition, the 
class encapsulates the precomputation of quantities that do not depend on the displacement vector, and it can return 
a pointer to a table of precomputations that do depend on a fixed displacement. For instance, for efficient evaluation 
of the translation of multipole expansions, Eq. (£7|), the class contains the (p + l) 4 values 

p-l Al A3-1 

Cg = ' (37) 

and an indirect-indexing array that specifies the sequence of 0^,i after the loops have been unrolled. In addition, 
a pointer to the (p + l) 2 precomputed values of Q k Y^ l {6, ip) can also be returned. Thus the transformation can be 
evaluated efficiently and with minimal overhead from loop-control variables. 

The computer memory requirements to store the precomputations that depend on spatial relationships can be 
greatly reduced for implementations that assume a regular lattice for spatial decomposition. For each level of the 
hierarchy, there will be a fixed number of displacements between the cells in that level and with cells on the parent 
level and the child level. The class sh_expansion encapsulates this efficiency at each request for a precomputation 
by searching a linked list of previous results and creating a new result only when no previous result exists. A 
similar scheme is also necessary for the compilation of multipole expansion coefficients for the lowest level of the 
hierarchy, Eq. (E^). Since the precomputations represent a significant amount of the memory requirements of the 
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overall simulation, these memory savings can greatly increase the number of lattice points used to represent the model 
magnet. 

The partitioning of space is accomplished through the class hbox. Instances of this class contain geometric infor- 
mation about the decomposition cell, expansions for both the potential and charges within the cell, as well as links 
to other cells. These links point to the encompassing cell, the encompassed cells, a list of same-level cells in the 
far field (the "interacting cells"), and a list of nearest-neighbor cells. The hboxs for each level of the decomposition 
hierarchy are held in a container class hlevel, and the entire hierarchy is maintained as a linked list. Our imple- 
mentation of hlevel is particular to the cubic-lattice decomposition of space, and irregular geometries have to be 
padded with empty space. (This is also true for our implementation of the classes VectorField and ScalarField 
from VectorField. h used throughout our numerical integration of the LLG equation.) 

APPENDIX 2 

A useful way to represent the thermal Landau-Lifshitz-Gilbert equation is in a form with the deterministic and 
stochastic parts separated 

dM(r < ,t) = B(M(ri,t))H do t(ri,i)dt + ViB(M(r i ,t))dW(r i ,t) , (38) 

where Hdot (r^ , t) is the deterministic part of the local field at and H„ (r^) = ^/edW is the stochastic part. 

The matrix B is given by 

/ a (Ml + Ml) —M z — aM x M y M y - aM x M z \ 
B(M(r,,t)) = — — M z — aM x M y a(Mj + M z 2 ) -M x ~aM y M z \, (39) 
+ a V -M y - aM x M z M x - aM y M z a (Af| + Af*) J 

where x, y, and z represent the Cartesian coordinates and the space and the tpae-dependence of the M^ have been 
omitted for clarity. The stochastic nature of the field results from the Wienerc-3~E3 process W (r^ , t) which has the 
properties 

(W„ (n, t)) = (W^ (r h t) WV ( r i> O) = (* ~ 0<W'<*M' ■ ( 40 ) 

The stochastic differential equation (^) can be treated numerically using first-order Euler integration. For small At, 
the deterministic integral is 

Idet(ri,i) =B(M(r i ,t))H det (r i> *)At . (41) 

The integral of the stochastic part, 

/t+At 
B(M(r 4 ,t'))dW(r 4 ,0 , (42) 

takes more consideration since it involves the product of the magnetization with the Wiener process. In such cases pi 
multiplicative noise, different methods for evaluating Eq. (M3) correspond to different Fokker-Planck equations .EaEZl 
There are an infinite number of ways to interpret Eq. (HjJ, but usually only the two extreme cases, the Ito and 
Stratonovich interpretations, are considered. The Fokker-Planck equation considered by Browncl only has the proper 
equilibrium properties when interpreted in the Stratonovich sense.0 This is complicated by the fact that numerical 
implementation of the stochastic integral is particularly convenient in the Ito interpretation. Then the discretized 
integral isE3 

I rt o(ri,*)=VeAtB(M(r i ,t))g(r i> *) , (43) 

where each component of g is a random number from a Gaussian distribution with zero mean and unit variance. 
Fortunately, changing interpretations can be accomplished through a correction discussed below. Then, Eq. (|43|) can 
be combined with Eq. ((4l|), and the result rearranged, to give the "impulse" during the integration step, Eq. (|ll|). 

Changing from the Stratonovich interpretation of the stochastic integral to that of Ito requires the addition of a 
deterministic term. Specifically, a Stratonovich interpretation of a multivariate-Langevin equation of the form 

dx = a (x, t) dt + b (x, t) dW , (44) 
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is equivalent to the Langevin equationS 

dx = (a (x, t)+c (x, t) ) dt + b (x, t) dW (45) 
in the Ito intcrpretationBE] where the components of the new drift term are 



c 



5 !>..%- <«> 



2 ^— ' ai/ 

For the present system, the additional drift term is readily found to be 

c(vi,t)= 6 M(r u t) , (47) 
(1 + a" 1 ) 

which is equivalent to the result given in Ref. |[ This drift term is always directed along the local magnetization 
density. The process of normalizing the magnitude of the spins during each integration step, Eq. (|l0|), is also directed 
along the magnetization, and essentially takes the correction into account. 
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3.6 nm 


M s 


1700 emu/cm 3 


Ar 


1.5 nm 


At 


5.0 x 10~ 5 ns 


7o 


1.76 x 10 7 Hz/Oe 


a 


0.1 



Table 1 . Parameters used in the micromagnetic simulations described here. 
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a) b) 

1.00 ns 1.05 1.10 1.15 1.20 ns 0.90 ns 0.35 1.00 1j05 1.10 ns 




FIG. 1. Magnetization reversal in the nanoscale magnet after a rapid reversal of the field. The magnets have a square 
cross section, but are shown in a one-quarter cut-away view. Both reversals are for H = 1800 Oe and T — 20 K. Here 
the z-component of the magnetization is shown, with red representing the metastable orientation and blue representing the 
equilibrium orientation, (a) Nucleation of both end caps, but at different times, (b) Nucleation of one single end cap that 
grows to reverse the entire magnet. 



a) b) 

0Oe -2325 -3235 -3300 -4000 OOe -2325 -3235 -3300 -4000 




FIG. 2. Magnetization reversal during a hysteresis loop of 1 ns at T — K. (a) The z-component of the magnetization with 
the same color scale as in Fig. 1 (b) The magnitude of the curl, |C|, with the sign taken from the value of C z . Red represents 
positive curl, blue negative curl, and green zero curl. The fields shown correspond to times 0.25, 0.375, 0.405, 0.45, and 0.5 ns, 
respectively. 
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FIG. 3. Hysteresis loops for the individual rectangular pillars at zero temperature and loop periods of 1, 2, and 4 ns. The 
change in the shape of the loops as the frequency is lowered indicates that the magnetization is not following the applied field in 
a quasi-static manner. The large tick marks on the upper horizontal axis indicate the times for the images of the magnetization 
and its curl, shown in Fig. H 
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FIG. 4. Energy density vs time for the single hysteresis-loop simulation at T = K with period 4 ns from Fig. ^. Assuming 
the energy is near its metastable minimum before the end caps start to propagate, the coercive field can be estimated from the 
local maximum in the total energy. 
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FIG. 5. Total energy density as a function of applied field for fields swept linearly with rate R at zero temperature. The 
estimation of the zero-rate coercive field H c estimated from the energy maximum is shown in the inset. The estimated value 
clearly depends on the rate of change of the applied field, but using linear fitting the static coercive field is found to be 
1979 ± 14 Oe. The error bars are estimated using the second derivative at the maximum. 
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FIG. 6. Mean switching time as a function of (a) applied field at T — 20 K, and (b) temperature at H = 1850 Oe. The 
switching time increases as (a) the free-energy barrier separating the metastable and stable orientations grows as the applied 
field is decreased, or (b) the thermal energy available for crossing the barrier decreases as temperature is lowered. 
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0.2 0.4 0.6 0.8 1.0 1.2 0.0 10.0 20.0 30.0 40.0 50.0 

t (ns) t (ns) 

FIG. 7. Probability of not switching, P no t(t), for (a) 200 switches of nanoscale magnetic pillars at H = 1800 Oe and T — 20 K, 
and (b) 2000 switches in a simple one-dimensional model nanopillars at H — 1000 Oe and T — 20 K. The theoretical forms, 
Eqs. @ and @, were fitted by matching their first and second moments to those of the simulation data. 




FIG. 8. Integration of weights to get P no t(t), Eq. (^). The solid lines separate regions where one (or two) nucleation events 
occur before switching, while the dashed curves are integration paths for constant switching time, t sv /to = 3/2 and 3. See the 
text for a full explanation. 
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FIG. 9. Schematic of two pillars, one with a growing region of magnetization in the equilibrium orientation at each end. 
The distance between pillars is d, their height is L, and the length of the reversed regions along the long axis are It and h, 
respectively. 




0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 



t (ns) t (ns) 

FIG. 10. (a) Probability of not switching, P no t(t), for nanoscale magnetic pillars in a four-pillar array at H = 1800 Oe and 
T — 20 K. The two pillars on the outside, and the two on the inside, are equivalent by symmetry. The isolated pillar data are 
the same as in Fig. m (b) Interactions between pillars as seen in the difference between P no t(t) for inside pillars with one and 
both neighboring pillars already switched. Here t=0 corresponds to the time of the most recent switch of a neighboring pillar. 
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FIG. 11. Schematic of the relationship between the lattice and the dual lattices, projected along the z-axis. The magnetization 
density M(rj) and the local magnetic field H(r-i) are known at the lattice sites of the simple cubic lattice (solid circles), the 
magnetic charge density pm is known at the dual lattice sites (open circles) which are located at the body center positions, and 
the magnetic surface charge density is known at the centers of squares defined by the surface lattice sites (open triangles). 



FIG. 12. Schematic of the hierarchical decomposition of space chosen for our implementation of the Fast Multipole Method. 
This eight-fold decomposition at each level is quite efficient given the underlying cubic lattice. 
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